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In  the  computation  of  discontinuous  solutions  of  hyperbolic  conservation  laws, 
TVD  (total-variation-diminishing),  TVB  (total- variation- bounded)  and  the  recently- 
developed  ENO  (essentially  non-oscillatory)  schemes  have  proven  to  be  very  useful. 
In  this  paper  two  improvements  are  discussed:  a  simple  TV  D  Runge-Kutta  type 
time  discretization,  and  an  ENO  construction  procedure  based  on  fluxes  rather 
than  on  cell  averages.  These  improvements  simplify  the  schemes  considerably 
especially  for  multi-dimensional  problems  or  problems  with  forcing  terms  Prelimi¬ 
nary  numerical  results  are  also  given. 


'Research  for  the  second  author  was  supported  by  the  National  Aeronautics  and 
Space  Administration  under  NASA  Contract  No.  NAS1-I8107  while  he  was  in 
residence  at  the  Institute  for  Computer  Applications  in  Science  and  Engineering 
(ICASE),  NASA  Langley  Research  Center,  Hampton.  YA  23065.  Also,  he  was 
supported  by  NSF  Crant  No  DMS85-03294,  ARO  ( J rant  No.  DA  AG29-85-K-0190, 
DARPA  Grant  in  the  ACMP  Program.  ONR  Grant  N(X)()1  I-86-K-069I .  and  N  ASA 
Langley  Grant  NAG  1-270 

l  lain  doebmoat  kf  b— ■  apiprmrwf 

lor  p*Kg*  Milan  <ad  adm  Ml 

dlstrfbuiftaa  to  idtaMgi  'f* 


□  □ 


I.  Introduction.  In  this  paper  we  are  interested  in  solving  the  system  of 
hyperbolic  conservation  laws 
d 

(1.1a)  «t  -f  ^2  /.  («)*,  =0  (or  =  g(u,  x,  t),  a  forcing  term) 

i=l 

(1.1b)  ti(x,0)  =  u0(x) 

Here  u  =  («j, . . .  ,um)T,  x  —  (xj,  xa, . . .  ,  x<j),  and  any  real  combination  of  the  Jaco- 

d 

bian  matrices  ^  has  m  real  eigenvalues  and  a  complete  set  of  eigenvectors. 

i'=i 

On  a  computational  grid  xy  =  j  •  Ax,  tn  =  nAl,  we  use  u"  to  denote  the 
computed  approximation  to  the  exact  solution  u(xy,  tn)  of  (1.1). 

We  also  use  the  abstract  form 

(1.2)  ut  =  C(u) 

in  place  of  (1.1a).  Here  £  is  a  spatial  operator. 

As  is  well  known,  the  solution  to  (1.1)  may  develop  discontinuities  (shocks,  con¬ 
tact  discontinuities,  etc.)  even  if  the  initial  condition  uo[x)  in  (1.1b)  is  a  smooth 
function.  Traditional  finite  difference  methods,  even  if  linearly  stable,  often  give 
poor  results  in  the  presence  of  shocks  and  other  discontinuities.  Recently  there 
has  been  a  lot  of  activity  geared  towards  constructing  efficient  finite  difference 
approximations  to  (1.1).  These  include  TVD  (total-variation-diminishing),  TVB 
(total- variation-bounded)  and  ENO  (essentially  non-oscillatory)  methods.  See,  e  g. 
[2J,  [3J,  (4),  (5),  [6],  [9],  (10J,  [12],  [13],  [14],  and  the  references  listed  therein.  Many 
of  the  ideas  can  be  traced  back  to  Van  Leer’s  work  in  [15],  [16). 

Usually,  rigoious  analysis  (e.g.  total- variation  stability,  convergence)  is  only 
done  for  the  scalar,  one-dimensional  nonlinear  case  (i.e.  d  —  m  =  1  in  (1.1)).  Some 
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partial  theory  (e.g.  convergence  for  first  order  monotone  schemes,  and  maximum 
norm  stability  for  higher  order  TVD  schemes)  exists  for  scalar  multi-dimensional 
problems  (d  >  1  in  (1.1)),  but  a  full  convergence  theory  for  multidimensional  non¬ 
linear  systems  appears  to  be  extremely  difficult.  However,  numerical  experiments 
for  multi-dimensional  problems  and/or  for  systems  of  equations,  using  direct  gener¬ 
alizations  of  TVD,  TVB  and  ENO  schemes  give  very  good  results.  Again,  see,  e.g., 
[2],  (4),  [6],  [10],  [11].  We  now  shall  confine  our  discussion  at  first  to  this  one  space 
dimension,  scalar  case.  Systems  and  multi-dimensional  problems  are  discussed  at 
the  end  of  Section  3. 


We  shall  always  use  conservative  schemes  of  the  form 


»7+‘=«j-a(4+j  -4_t),  a  = 


with  a  consistent  numerical  flux 


fj+i  =  y+fc);  /(« - «)  =  /(«) 


in  order  to  guarantee  that  any  convergent  bounded  a.e.  subsequence  has  as  its  limit 
a  weak  solution  of  (1.1),  (Lax-Wendroff  Theorem  [8]),  i.e.  we  construct  so-called 
“shock  capturing  methods”. 


The  total  variation  of  a  discrete  scalar  solution  is  usually  defined  bv 


TV(u)  =  /  Juy+i  -  Uy| 


We  say  the  scheme  is  TVD  if 


TV(un+l)  <  TV(un) 


for  some  fixed  B  depending  only  on  TV(u°),  and  for  all  n  and  A t  such  that  0  < 
nAt  <  T. 

A  nice  theoretical  advantage  of  all  TVD  or  TVB  schemes  is  that  they  have  con¬ 
vergent  subsequences  as  Ax  — ►  0,  and,  if  a  further  “entropy  condition”  is  satisfied, 
then  they  are  convergent.  See,  e.g.  [3] . 

The  formal  “order  of  accuracy”  in  this  paper  is  in  the  sense  of  local  truncation 
errors,  i.e.  if  local  truncation  error  is  0(Az’’+l)  in  smooth  regions,  we  say  the 
scheme  is  (formally)  r-th  order  accurate.  See,  e.g.  [2]. 

There  are  many  TVD  schemes  constructed  in  the  literature  (e.g.  [2],  [3),  [9], 
[10],  [14]).  In  [10],  TVD  schemes  of  very  high  spatial  order  (up  to  15th  order) 
were  constructed.  These  schemes  can  be  used  for  steady  state  calculations  (e.g. 
implemented  with  the  TVD  Runge-Kutta  type  time  discretizations  with  large  CFL 
numbers  in  [13])  or  for  time  dependent  problems,  equipped  with  a  multi-level  TVD 
high  order  time  discretization  in  [13]  or  with  a  Runge-Kutta  type  TVD  high  order 
time  discretization  in  Section  2  of  this  paper.  These  are  perhaps  the  highest  order 
TVD  schemes  existing  at  present.  However  the  definition  of  total  variation  (1.5) 
implies  that  these  methods  must  degenerate  to  first  order  accuracy  at  extremacy.  A 
TVB  modification  of  such  schemes  which  recovers  global  high  order  accuracy  even 
at  critical  points  is  obtained  in  [12]. 

The  above  mentioned  TVD  and  TVB  schemes  use  a  fixed .  wide  stencil  (for  the 
15th  order  scheme,  the  stencil  is  17  points  wide),  thus  restricting  the  advantage  of 


going  to  higher  order  through  smearing  of  discontinuities  and  resulting  degradation 
of  the  accuracy.  Numerically  we  observed  that  third  order  schemes  work  quite  well 
[12],  but  we  lost  accuracy  in  a  fairly  large  region  near  discontinuities  by  using  a  fifth 
order  method.  Recently  Harten,  Osher,  Engquist  and  Chakravarthy  constructed 
ENO  schemes  which  are  of  globally  high  order  accuracy  in  smooth  regions  and 
which  use  adaptive  stencils,  thus  obtaining  information  from  regions  of  smoothness 
if  discontinuities  are  present.  These  methods  achieve  high  order  accuracy  right  up 
to  discontinuities.  Analysis  and  numerical  experiments  are  found  in  [5],  [6],  [4]. 
At  present,  a  convergence  theory  (e.g.  TV  boundedness)  for  ENO  schemes  is  still 
unavailable. 

There  are  two  natural  directions  in  which  to  simplify  the  ENO  or  TVD,  TVB 
schemes,  especially  for  multi-dimensional  problems  or  problems  with  forcing  terms: 

(1)  Time  discretization.  Usually,  semi-discrete  (method  of  lines)  versions  of 
ENO  or  TVD,  TVB  schemes  are  much  simpler  than  the  fully  discrete  ones.  There 
are  then  mainly  two  ways  to  discretize  in  time.  One  is  of  Lax-Wendroff  type,  i.e., 
by  using  ut  =  -/*,  u«  =  (/'/*)*, •  •  ■ ,  «y  +  l  =  «y  +  A *(««)*+  ^(«m)* + 
and  then  by  discretizing  the  spatial  derivatives.  Many  second  order  TVD  schemes 
(e.g.  Harten’s  in  [2]),  and  the  ENO  schemes  in  [5],  [6],  [4],  used  this  type  of 
time  discretization.  The  main  disadvantages  to  the  procedure  is  that  it  is  com¬ 
plicated  to  program,  especially  for  multi-dimensional  problems  with  forcing  terms. 
One  can  see  this  by  writing  out  a  third  order  approximation  to  the  equation 
ut  4-  /,(«)„  +  =  g(u,ii,x3,/).  Moreover  it  is  not  easy  to  prove  that 

this  results  in  a  TVD  or  TVB  method,  even  if  the  original  method-of-Iines  ODE 
and  its  Euler  forward  version  are  both  TVD  or  TVB.  The  numerical  results  have 
proven  satisfactory,  but,  speaking  theoretically,  only  second  order  in  time  TVD  or 


TVB  schemes  exist  of  this  type.  Another  way  to  discretize  in  time  is  to  use  a  multi¬ 
level  or  Runge-Kutta  type  ODE  solver.  This  is  much  simpler  to  program  than  the 
Lax-Wendroff  type  of  discretization  for  multi-dimensional  problems  or  for  problems 
with  forcing  terms,  so  it  is  widely  used  for  numerically  implementing  a  method  of 
lines  approximation.  However,  usually  only  linear  stability  analysis  is  available  in 
the  literature,  which  is  certainly  not  enough  for  our  purpose  since  linear  stabil¬ 
ity  does  not  imply  convergence  if  shocks  or  other  discontinuities  are  present.  This 
is  particularly  true  for  ENO  schemes  which  use  moving  stencils.  Linear  stability 
analysis  is  based  on  the  fact  that  the  stencil  is  fixed  and  the  error  accumulates  in 
a  predictable  pattern,  hence  it  does  not  apply  to  ENO  schemes  at  all.  For  these 
reasons  we  consider  TVD  time  discretizations  .  In  [13],  a  class  of  multi-level  TVD 
time  discretizations  were  constructed  and  analysed,  (numerical  results  can  be  found 
in  [12]).  However,  for  easy  starting  and  for  storage  considerations,  one  step  Runge- 
Kutta  type  schemes  are  preferable  to  multi-level  methods.  In  Section  II  of  this 
paper  we  present  a  class  of  high  order  TVD  Runge-Kutta  type  time  discretizations. 

(2).  Avoiding  the  using  of  cell-averages.  The  ENO  schemes  constructed  in  [5], 
[6],  [4]  are  for  cell-averages  but  involve  point  values  as  well.  Hence  a  reconstruction 
procedure  is  needed  to  recover  point  values  from  cell  averages  to  the  correct  order, 
which  can  be  rather  complicated,  especially  in  multi-dimensional  problems.  It  is 
desirable  to  use  the  moving-stencil  idea  directly  on  fluxes  to  get  ENO  schemes 
without  using  cell-averages.  In  Section  III  of  this  paper  a  class  of  such  ENO  schemes 
is  constructed. 

Some  encouraging  numerical  results  obtained  by  using  schemes  constructed  in 
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this  paper  are  included  in  Section  IV. 
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We  conclude  these  introductory  remarks  by  noting  that  R.  Sanders  [17]  has 
recently  devised  third  order  accurate  TVD  methods  which  degenerate  to  second 
order  at  extrema.  He  defines  the  variation  of  the  numerical  solution  as  the  variation 
of  an  appropriately  chosen  piecewise  parabolic  interpolant.  The  numerical  results 
are  very  good.  However  this  technique  has  no  method  of  lines  analogue,  so  we  omit 
it  from  our  present  discussion. 


II.  High  Order  Runge-Kutta  Type  TVD  Time  Discretizations.  Define 


w  =  T(u)  =  (/  +  Af£)(u) 


where  T  and  L  are  nonlinear  discrete  operators,  £  is  a  r-th  order  discrete  approxi¬ 
mate  to  the  spatial  operator  £.  in  (1.2): 


L(n)  =  £(«)+  0(A*r) 


if  u  is  smooth. 


Our  goal  is  to  get  a  fully  r-th  order  approximation  to  the  differential  equation 
(1.2)  of  the  form 


u"+l  =  S(un) 


(The  operator  S  depends  on  T).  This  means  that  if  u(x,<)  is  an  exact  smooth 
solution  of  (1.2),  then 


u(xy,r+1)  -  S(un)y  =  0(Axr+1 ). 


We  also  want  the  scheme  to  be  TVD: 


TV(S(u))  <  TV(T(u )) 


under  suitable  restrictions  on  At  (or,  equivalently,  on  the  CFL  number  A). 


vOT 
:>.*/- 


I 


v.v.v 


We  call  a  time  discretization  (2.3)  r-th  order  TVD  if  it  satisfies  (2.4)  and  (2.5). 
If  the  spatial  operator  T  in  (2.1)  is  TVD  or  TVB: 

(2.6a)  iV(T(«))  <  TV(u) 

or 

(2.6b)  TV(T(u))  <  TV(u)  +  MAt 

for  0  <  M  uniformly  bounded  as  At  — +  0,  then  the  fully  discrete  high  order  scheme 
(2.3)  is  TVD  or  TVB,  owing  to  (2.5). 

In  [13],  a  class  of  multi-level  type  high  order  TVD  time  discretizations  was 
constructed.  Numerical  experiments  in  [12]  were  very  promising  .  But  there  are  two 
disadvantages  of  multi-level  type  methods:  (i)  for  an  m-th  level  method  the  first  m  — 
1  levels  have  to  be  calculated  by  other  methods  to  the  same  order  of  accuracy  (e.g. 
by  using  Taylor  series  expansions);  (ii)  we  have  to  store  all  m  level  datas,  creating  a 
rather  large  storage  requirement,  stretching  up  to  and  beyond  the  limits  of  present 
day  computers  for  physical  problems  arising  e.g.  in  computational  aeronautics.  At 
present,  Runge-Kutta  type  methods  are  more  often  used  in  discretizing  the  method- 
of-lines  than  are  multi-level  methods.  Since  the  former  consists  of  one-level  methods, 
they  are  self-starting  and  reduce  storage  requirements  significantly.  In  the  following 
we  will  analyze  the  nonlinear  stability  (TVD)  of  a  class  of  such  methods. 

Assume  (2.1)  is  TVD  (or  TVB)  under  a  suitable  CFL  restriction 
(2.7)  A  <  A0 

We  may  also  need  an  approximation  to  -L  to  the  spatial  operator  which 
we  take  to  satisfy 


(2.8) 


w  =  T(u )  =  (/—  AtL)(u) 


where  (2.8)  is  TVD  (or  TVB)  under  the  same  CFL  restriction  (2.7). 


As  an  example,  a  very  simple  first  order  (r  =  1)  TVD  approximation  to 


ut  =  ux  =  £(u) 


is  obtained  via  simple  upwind  differencing: 

nu\  _  +  Ax)  -  tt(ar) 

Ax 

This  scheme  (2.1)  is  TVD  if  A  =  ££  <  A0  <  1. 


The  approximation  L{u)  is  defined  as 

r/..x  u(x)  —  u(x  —  Ax) 

LW)  ~ - Al - 

and  (2.8)  is  TVD  again  for  A  =  <  A0  <  1. 


This  procedure  (2.8),  (2.9)  easily  generalizes  to  any  conventional  TVD,  TVB, 
or  ENO-approximation  (2.1)  satisfying  (2.2). 


The  general  explicit  Runge-Kutta  method  (we  use  explicit  methods  to  avoid 
solving  nonlinear  equations)  for  (2.1)  is 


(2.10a) 

(2.10b) 


«(*)  =«(0)  + At  £c, •*£(«<*>),  *  =  1,2,..., 

k=0 

u<°>  =  u",  u(m)  =  un+l 


If  the  operator  L  also  depends  explicitly  on  t,  as  is  the  case  when  the  forcing 
term  g  in  (1.1a)  depends  explicitly  on  t,  or  when  we  have  time-dependent  boundary 
conditions,  the  general  explicit  Runge-Kutta  method  takes  a  more  complicated  form 


(2.11a) 


where 


u(»)  _  yfo)  +  ^  c,fcL(u*fc),  f(0)  +  dfcA/) 
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For  the  classical  4-th  order  Runge-Kutta  methods,  the  constants  c,*  in  (2.10) 
are  all  non-negative.  However,  /?,*  in  (2.12)  may  well  be  negative.  In  order  to  obtain 
TVD  we  apply  a  trick  used  in  [13],  i.e.  we  replace  L  in  (2.12)  by  L  in  (2.8)-(2.9) 
whenever  /?,•*  is  negative.  Now  (2.12)  becomes  a  convex  combination  of  TVD  (or 
TVB)  operators  under  the  CFL  restriction 


A  <  A0min-rp^- 


(2.13) 


and  we  easily  get  the  following 


PROPOSITION  2.1.  Scheme  (2.12)  is  TVD  under  the  CFL  restriction  (2.13), 
if  L  is  replaced  by  L  when  is  negative. 

REMARK  2.1.  The  previous  proposition  may  be  put  into  a  more  general 
framework  as  follows.  The  TV  in  (2.5)  and  (2.6) (a, b)  may  be  replaced  by  (?(«), 
any  convex  mapping  into  the  non-negative  real  line,  where  u,  S(u),  T{u )  belong 
to  a  Banach  space  of  functions  BAl.  Also  if  u  is  a  smooth  solution  of  1.2,  then 
u(x,tn+1)  —  S(u(x,tn)}  =  0((Aj)r+1).  The  statement  in  the  proposition  can  then 
be  replaced  by:  G(um)  <  G(u°);  moreover  the  formal  order  of  accuracy  is  still  r. 

Now  our  goal  is  to  choose  the  a ,•*  and  /?,•*  such  that  (2.12)  is  of  the  highest 
possible  order  and  such  that  the  CFL  restriction  (2.13)  is  optimal.  We  would  also 
like  to  minimize  the  number  of  negative  /?,*’  s  in  order  to  reduce  the  computational 
work  involving  L. 

One  easy  way  to  do  this  is  to  use  ?  standard  Runge-Kutta  method  and  then 
rewrite  it  in  the  form  (2.12)  to  get  a,*  and  /?,*.  Unfortunately,  most  classical 
Runge-Kutta  methods  lead  to  small  CFL  numbers  in  (2.13)  as  well  as  negative 
Jilt's.  Hence  the  best  way  is  to  consider  (2.12)  directly.  Straightforward  but  tedious 
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Taylor  expansions  and  an  analysis  of  possible  parameters  (which  we  omit)  leads  us 
to  the  following  results: 


i)  Second  order  case,  m  =  1. 


For  accuracy 


(2.14) 


o; 20  —  1  —  O21 
@20  ~  1  ~  2^Io  ~  a21.^10 
=  oh 


Pio,  <>2i  are  free  parameters. 


It  can  be  verified  that  the  “optimal”  scheme  (considering  CFL  restriction  (2.13) 


and  whether  L  appears)  is 


(2.15) 


1 1  =  u(0*  +  A  tL(u^) 
u(2)  =  l«(o)  +  L^tL{u^)  . 

CFL #  =  1 


Notice  that  L  does  not  appear  in  (2.15).  This  is  equivalent  to 
(«<*>  =  u<°>  +  A<L(u<°>) 

(216)  l  u*2>  =  u<°>  +  |A tL(uW)  +  1A tL(u^) 


which  is  the  classical  Heun’s  method  or  modified  Euler  method  111. 


ii)  Third  order  case,  m  =  2. 


For  accuracy 


p32  — 


6P[0io-P) 


(2.17) 


.3  _  *  —  <*33  010  031—  P  033 

J3l  _  - - - 

>ho  =  1  —  «3l  PlO  -  q32  P  -  thl  —  ‘^32 

@20  —  P  ~  <>21  PlO  ~  P‘21 


(>2i -  q30,  031,  fiio  and  P  =  @2o  +  021  dio  +  $21  are  free  parameters.  The  solution 


(2.17)  is  written  in  convenient  inductive  form. 


*|i  >«|t  ii.  4U  iH  *  »*« 


Extensive  searching  leads  to  the  following  preferred  scheme: 


W*\v 

v'v 

§83$ 


(2.18) 


u(i)  =tt(°)+A/I(tt(°>) 

ul3)  =  |u(°)  4-  4-  ^AtL(u^) 

„(3)  =  1„(0)  +  |u(3)  +  |AtI(«(2)) 
CFI#  =  1 


Hast? 


Notice  that  £  does  not  appear  in  (2.18).  Also  in  computing  u **)  we  only  need 
L(u(*~1)),  so  there  is  no  need  to  store  the  previous  L(u^),  lowering  the  storage 
requirement  significantly. 

(2.18)  is  equivalent  to 


-I-.1 

.’wv^vy.'j 

• 

>->>v 

•v  a 


(2.19) 


„(i)  =  «(°)  +  AtL(uW) 

«(*>  =  «(°)  4-  ±Af£(u<°>)  4-  \AtL(u^) 

«(3)  =  «<°)  +  idf£(«<°))  +  iAf£(u<^)  +  §Af/;(u<3>) 


JU*  -*  F-  A  ft. 

V’  rVl 


We  have  been  unable  to  identify  (2.19)  with  any  of  the  “classical”  third  order 
Runge-Kutta  methods.  On  the  other  hand,  the  classical  third  order  Runge-Kutta 
methods  in  [7j,  when  written  in  equivalent  forms  (2.12),  lead  to  negative  /?,■*  and 
small  CFL  numbers  (2.13),  and  are  hence  inferior  to  (2.19). 


(iii)  Fourth  order  case,  m  =  3 


For  accuracy  we  get  a  system  of  7  equations  with  16  unknowns,  so  there  are 
9  free  parameters.  Unfortunately  this  time  the  solution  is  not  easily  obtained  in  a 
convenient  form.  In  [1]  a  general  solution  with  two  parameters  is  given  for  the  form 


/v 


.vr-r^\ 

N  >  .  ' 

Xi/V 


^7" 

V  /  V 

V  v  ■. 

>,v.v. 

•  •/». 
V  /  V 


(2.10).  We  can  certainly  rewrite  it  in  the  form  (2.12).  Extensive  searching  seems  to 
indicate  that  we  cannot  avoid  negative  constants  /?,*  this  time.  The  classical  fourth 


S*S*!*vS 

,\-VA 

•'•V/ 


order  Runge-Kutta  method  can  be  written  in  the  form  (2.12)  as 


(2.20) 


■  tiO)  =  «(°)  +  IA  <I(u<°>) 

«(3)  =  l„(o)  -  1A tZ(«<°>)  +  +  §A tL(u^) 

.  u(3)  =  Iu(o)  _  lAfZ(ti<°>)  +  |ti^)  -  ±A tL{uW)  +  §u(2)  +  AfZ(u(2>) 
«(«>  =  |«(‘)  +  iAt^u*1))  +  ±«(2>  +  +  iA  tL(uW) 

. CFL#  = | 


Notice  that  we  have  to  compute  Z(«(°^)  and  Z(u(1)).  If  we  use  the  more 
awkward  definition  of  « ^ 


(2.21) 


6431  /0>  18769  *  (0».  18769  /ji  137  ~  rj>  137  /a>  .  o » » 

«(3)  =  rrrrr  «<°>-7rrrrrAf  i(«(0))+^^«(1)-— At  L(«(1))+— u(2,+A<L(u<3>) 


80000 


160000 


80000 


400 


200 


then  the  CFL#  can  be  raised  slightly  to  *§*. 


Another  fourth  order  scheme  with  a  slightly  larger  CFL  #  is: 

’  ul1)  =  «(°)  +  §A<L(ti<°>) 

=  |«(°)  -  £a<Z(«<°>)  +  |u<l>  +  |AtL(«<1)) 

,2  22.  «(3)  =  5&«(o)  -  »^(«(o)) + s«(i) "  S8^(«(1)>+ 

('  ]  +^U(2)  +  |A/L(«(2)) 

„(4)  _  l„(o)  +  lu(i)  +  1A tL(u<l>)  +  £u<3>  +  £A/L(«<3>) 
k  CFL #  =  0.87 

We  still  need  to  compute  Z(u*0*)  and  Z(u^*). 


(iv)  Fifth  order  case,  m  =  5. 


We  simply  write  out  the  form  (2.12)  corresponding  to  a  fifth  order  method 


given  on  page  143  of  [7]: 


(2.23) 

u*1)  =  +  |Af£(u*°)) 

UW  =  |«(°)  +  |u^  +  f  A  tL(u^1^) 

u(3)  _  |u(o)  _  LAll(u(o))  +  Iu(i)  _  JLA t'L(u^)  +  ±«(2>  +  ±A tL(uW) 
u<4)  =  l„(o)  _  A  A  t'L(uW)  +  Ail*1)  -  }}  At'L(u^)  +  +  f  AfL(u<2>)  + 


u 


(5)  _89S37_  (0)  ,  2270219  Afjr 

1  2880000“  “  40320000^^ 

-1£LL„(2)  .  miAfI(U(2>) 
12000“  ^  2800^llJ'“  ' 


T  200 


«(•)  =  |ti(°)  +  -fruM  -  £A tL(u^)  +  +  | A tL(uW)  +  «<5>  +  ^AtL(n^) 

.CFL#=^ 


Notice  that  we  need  to  compute  Z(td0)),  L(u * 1 ) ),  L(u^). 


III.  A  Simplified  Version  of  ENO  Schemes.  To  use  the  Runge-Kutta 
type  TVD  time  discretizations  in  Section  II,  we  must  have  a  spatial  discrete  opera¬ 
tor  (2.1)  to  start  with.  Theoretically  one  would  like  to  use  a  TVD  or  TVB  operator 
T  satisfying  (2.6),  because  then  the  full  scheme  (2.3)  would  be  TVD  or  TVB.  But 
as  indicated  in  section  I,  the  existing  high  order  TVD  or  TVB  schemes  may  smear 
discontinuities  and  pollute  the  solution  (i.e.  we  may  not  get  high  order  accuracy  in 
a  fairly  large  region  near  discontinuities),  due  to  the  fixed,  wide  stencil.  The  ENO 
schemes  constructed  in  [5],  [6],  [4]  are  very  promising  experimentally  and  appealing 
conceptually,  but  the  fact  that  they  use  cell-averages  as  well  as  point  values  via  a  re¬ 
construction  procedure,  and  that  they  were  implemented  using  a  Lax-VVendroff  type 
time  discretization,  makes  them  rather  complicated  to  program,  especially  in  multi¬ 
dimensional  problems,  or  problems  with  forcing  terms.  The  Runge-Kutta  type  TVD 
time  discretizations  in  Section  II  equipped  with  semi-discrete  ENO  schemes  will  sim¬ 
plify  them  in  many  cases  (although  there  is  no  rigorous  theory  concerning  TVB  of 
semi-discrete  ENO  schemes  or  their  Euler  forward  version,  analysis  in  many  cases 
and  numerical  experiments  strongly  support  that  the  total  variation  increase  at  each 


step  is  0( Axr)  for  a  r-th  order  ENO  scheme,  see  [6).  Hence  the  full  scheme  (2.3) 
in  this  case  should  also  be  TVB).  In  this  section  we  further  simplify  non-oscillatorv 
methods  by  deriving  a  version  of  ENO  schemes  using  only  fluxes,  not  cell-averages. 


We  start  with  a  simple  first  order  monotone  Lax- Friedrichs  type  of  scheme.  If 
we  define 


/+(°)  =  «(/(«)  +  au).  /  (u)  =  »(/(«)  -  au) 


where  a  >  max  (  is  a  constant,  then  clearly 


/+'(«)  >  o,  r»  <  0 

/+(«) +r  («)  =  /(«) 


The  Lax- Friedrichs  scheme  is  simply  (1.3)  with  the  numerical  flux  defined  by 


/**  =  //+/; 


Taylor  expansion  reveals  the  existence  of  constants  a2.  a*, . . . ,  a2rn_3 . such 

that  if 

f*~ 1  k 

(3.5)  =  /,t  j  +  £  +  0(Aram+l) 

fc=l 

then  the  scheme  (1.3)  will  be  2m-th  order  accurate  in  space  in  the  sense  of  (2.2) . 
For  example,  aa  =  - ± 

In  light  of  (3.4),  it  is  natural  to  require 

<3-«>  /,+!  =  /;+,+ /-+i 

and  to  define  the  positive  flux  f++i  and  the  negative  flux  f~+l  (in  the  meaning  of 
(3.2))  separately. 
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For  accuracy,  we  require  ft  x  and  /.  A  to  satisfy  (3.5)  separately: 

•  '  2  ^  '  1 

m— 1  o2  fc 

1  4*}  «  /,*}  +  E  “>^*“(a^E/±)y+*  +  <W"+1) 


We  achieve  (3.7)  by  using  polynomial  interpolants  p*  of  f±  to  the  correct 


order: 


Pj*(;r)  =/*(«(*)) +  0(A*3"*+1) 


near  x  =  xj+a,  then  define 


(3-9)  /£x  =Pf  (*,*+*)  +  E  fl^(Ax)2fc(^pJ)i=i 

fc=i  * 

Clearly  if  (3.8)  is  true,  then  the  fluxes  A  defined  by  (3.9)  will  satisfy  (3.7). 

It  is  in  constructing  the  interpolating  polynomials  pf(x)  that  we  use  the  ENO 
moving  stencil  idea:  in  order  to  achieve  (3.8)  pf(x)  can  be  polynomials  of  degree 
2m  interpolating  f±(u(x))  at  any  2m  +  1  points  near  Xj+i-  We  use  the  ENO  ideas 
in  [5],  [6]  and  [4]  to  choose  the  2m  +  1  points  automatically  from  the  smoothest 
possible  region,  but  start  with  the  correct  one  according  to  (3.4). 

The  algorithm  can  be  written  as  follows:  For  constructing  p+(x): 

(3.io)  (i)  C„  =  *!2!«  =  i.  Q{°'M  =  f+M 

(2)  Inductively,  assume  we  have  1 1 .  I^11  and  then  we  coni- 

pute  the  n-th  divided  differences  of  /+(u(x)): 


(3.11a) 

(311b) 


a(n)  =  /+[u(x  . »(x 

min  f  * 

fe<n'  =  /+[u(xfc,  „-!)_,) . «(xfc, 
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(b)  The  second  order  scheme  here  is  similar  to  the  usual  minmod  second  order 
TVD  scheme,  [9]  except  that  the  minmod  function  is  replaced  by  choosing  the  value 
closer  to  zero  (we  omit  the  details  of  the  derivation  here);  this  is  still  a  TVD  method 

(c)  For  the  linear  equation  ut  +  aux  =  0  (the  scheme  is  still  nonlinear!),  the 
schemes  here  are  equivalent  to  the  ENO  schemes  in  [6]  using  the  primitive  function 
reconstruction,  except  for  a  possible  difference  in  the  choice  of  stencil.  The  ENO 
schemes  in  [6]  choose  the  stencil  according  to  the  divided  difference  table  of  cell 
averages  of  u,  rather  than  that  of  /*  (u)  here.  We  again  omit  the  details  of  derivation 
here.  Since  the  ENO  schemes  in  [6]  worked  so  well  numerically,  and  our  simplified 
schemes  here  are  equivalent  to  those  in  [6]  for  linear  equations,  we  expect  ours  to 
work  as  well.  For  preliminary  numerical  results  see  section  IV. 

As  mentioned  before  there  is  at  present  no  rigorous  theory  about  TVB  of  this 
type  of  ENO  schemes.  Here  we  make  two  observations  for  our  schemes  along  these 
lines: 


(1)  For  smooth  solutions  all  the  divided  differences  (3  11)  should  be  bounded 
(by  the  maximum  norm  of  the  n-th  derivative  of  f±.  times  some  constant)  So  if 
we  use 

(3.16a)  c*n)  =  min(|c('l*|,  )  sign  (c*n| ) 

or 

(3.16b)  c(n)  =  min(|c*n||,  A/,n| Arn-3 )  sign  (c1'*1 ) 

in  the  place  of  c(n*  in  (3.14)  for  n  >  2.  where  ,\/(n|  are  constants  which  are  related 
to  the  maximum  norm  of  the  n-th  derivative  of  in  initial  smoot h  regions,  it  due*, 
not  affect  the  accuracy  in  regions  of  smoothness  We  then  get  a  TVB  scheme  We 
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can  easily  see  that  our  flux  with  (3.16)  satisfies 

(3»7a)  f,+  i=fJ?F  +  d>+i 

where 


wwrrwuw 


$Vv.\*v 

\w:/- 


(3.17b)  |cy+x|  <  A/A  I2 

with  the  constant  A/  depending  on  the  A/^n^,  and  /TVDa  is  the  second  order  TVD 

*  '  I 

flux  mentioned  in  Remark  (3.1b)  above.  Now  (3.17)  clearly  implies  TVB  of  the 
scheme.  Numerically  we  do  not  see  any  essential  difference  by  using  or  not  using 
(3.16),  hence  we  strongly  believe  that  ENO  schemes  without  (3.16)  are  also  TVB, 
at  least  for  most  practical  problems. 

(2)  The  approach  in  (3.1)  -  (3.15)  is  not  the  only  possible  one  which  can  be 
used  to  construct  an  ENO  scheme  based  on  interpolating  fluxes.  At  an  early  stage 
of  our  current  work  we  used  another  approach,  starting  from  f+  and  f~  in  (3.1), 
then  for  constructing  Py  (x): 

=  =  g(+l)u)  =  /+(^)  +  /+l«Uy-.)-«(^)](x-xy); 

(2)  k  (3),  same  as  the  procedures  (2)  and  (3)  above  in  (3.10), 

Similarly,  for  constructing  p ~  ( x ) : 

(1)  =  ;  +  1.  ^Ln(x)  =  /-(«,)  +  /_l«(xy).  (xJ  +  i)](x  -  x,); 

(2)  and  (3),  same  as  the  procedures  (2)  and  (3)  above  in  (3.15). 

Then  we  take  p}{i)  —  pf  (x)  +  p~ (x),  and  write  our  scheme  as 
(3.18)  ti?+l  =  «;  -  A/(^p,(x)),=If 
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Notice  that  the  scheme  (3.18)  is  simpler  than  (3.9)  -  (1.3).  The  only  trouble  is  that 
it  is  not  in  conservation  form  (1.3).  However,  if  we  use  (3.16),  then  it  is  easily  seen 
that  (3.18)  can  be  written  as 


(3.19) 


where  |cy|  <  A/,  and  /f^  is  the  first  order  Lax- Friedrichs  flux  (3.-1).  We  can  call 
such  schemes  “essentially  conservative”  because  the  most  important  property  of  a 
conservative  scbeme-the  conclusion  of  the  Lax-Wendroff  theorem  in  [8]  -  is  still 
valid.  Since  the  scheme  deviates  from  a  first  order  monotone  scheme  (not  only  a 
TVD  scheme)  by  M Ax3,  we  have  even  a  stronger  theory  than  before  -  we  have 
the  entropy  condition,  hence  full  convergence  (not  just  of  a  subsequence),  and  also 
convergence  in  multi-dimensional  scalar  problems  i.e.  we  have  every'  convergence 
property  first  order  monotone  schemes  have.  Unfortunately  numerical  experiments 
indicate  that  in  some  cases,  (3.18)  is  inferior  to  the  fully  conservative  (3.9)  -  (1.3). 
An  illustrative  example  is  to  compute  the  Riemann  problem  for  Burgers'  equation 
ut  +  uu,  =  0  with  a  moving  shock  (e  g.  u\en  =  Hrjght  =  — 1.)  using  the  fifth 
order  versions  of  (3.18)  and  (3.9)  -  (3.15),  (1.3),  equipped  with  a  fifth  order  multi¬ 
level  TVD  time  discretization  in  [13].  The  procedure  (3.9)  -  (3.15),  (1.3)  gives  good 
results,  with  or  without  (3.16);  while  without  (3.16).  the  non  conservative  (3.18) 
gives  the  wrong  shock  location.  With  (3  16)  the  shock  location  becomes  correct,  but 
the  mechanism  that  enforces  this  causes  a  rather  severe  smearing  of  the  shock.  For 
these  reasons  we  abandoned  the  simple  and  theoretically  pleasing  version  (3.18) 

Finally  let  us  point  out  that  the  Lax- Friedrichs  building  block  is  only  con¬ 
venient  one:  we  may  also  use  other  monotone  or  £- fluxes  (see.  e  g  [it)])  as  our 
building  blocks  Of  course  it  is  not  always  possible  to  assoicate  /+  and  f~  as  in 
(3  1)  with  each  ff-flux  such  that  (3.2).  (3  3).  (3  t)  is  valid,  but  careful  inspection 
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reveals  that  we  do  not  need  to  use  the  values  of  /+  and  /  -  only  their  divided 

differences.  For  each  E-flux  /iy+x  we  can  define 

(3.20)  4/y+ a  =  /y+i  ~  ^y+  4  >  ^/y+x  ~  hi+k  ~  fi 

where  d/t  A  and  df~  x  replace  the  first  (undivided)  differences  -  /,+  and 
—  f~ .  Hence  we  can  just  use  the  divided  difference  tables  of  df+  and  df~  in 
place  of  the  divided  difference  tables  of  /+  and  f~  in  constructing  py  and  p~,  and 
define  ff,  f~  in  any  consistent  way  such  that  /+  +  =  f, ,  e.g.  /t  =  /,,  /“  =  0. 

By  (3.20) 

(3.21)  df++i  +  d/~+ 1  =  /y+ 1  -  /y 

hence  if  p+  and  p"  use  the  same  stencil  then  pt(x)  +  p~ (x)  is  a  polynomial  inter¬ 
polating  /(u(x)),  thus  accuracy  is  guaranteed  with  (3.9)  -  (3.6);  if  the  stencils  are 
different,  say  p+  has  the  same  stencil  as  pj  but  p+  does  not,  then  it  is  easy  to  show 
that  pt  -pt  is  a  sum  of  r-th  order  undivided  differences  of  dff  A  (r  is  the  order  of 
the  polynomials  p+,  p~)  hence  as  long  as  these  are  0(Axr+1)  (valid  if  the  E-flux 
/ly+i  is  smooth  up  to  order  r)  we  still  have  the  correct  accuracy. 

The  reason  one  might  consider  general  E-fluxes  as  building  blocks  is  that  the 
Lax- Friedrichs  flux  is  considered  to  be  too  dissipative.  While  the  first  order  Lax- 
Friedrichs  scheme  is  much  inferior  to  upwind  schemes  (e  g.  to  Godunov's  or  the 
Engquist-Osher  schemes),  our  numerical  experiments  show  that  higher  order  ENO 
schemes  using  Lax-Friedrichs  building  blocks  work  quite  well  (although  they  are  still 
slightly  inferior  to  the  same  order  ENO  schemes  based  on  upwind  building  blocks, 
the  difference  is  much  smaller  than  that  in  the  first  order  case).  The  advantage  of 
the  Lax-Friedrichs  flux  is  that  it  is  C°° ,  hence  the  ENO  schemes  based  on  it  have 
full  high  order  accuracy.  On  the  other  hand,  most  other  E-fluxes  Godunovs, 
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Engquist-Osher’s,  Entropy  condition  satisfying  version  of  Roe’s,  etc.  -  are  not 
smooth  at  sonic  points  (points  at  which  f'{u)  =  0),  hence  ENO  schemes  based 
on  them  using  the  methodology'  of  this  paper  will  lose  accuracy  at  sonic  points. 
Although  we  may  overcome  this  by  smoothing  those  Effluxes  at  sonic  points,  in 
most  cases  the  simple  Lax- Friedrichs  building  block  should  be  good  enough. 


Problems  in  multi-dimensions  are  approximated  by  applying  the  procedure 
described  in  (3.1)  -  (3.15)  or  its  generalization  (3.20),  (3.21)  to  each  of  the  terms 
in  (1.1a).  The  Range-Kutta  methods  devised  in  section  2  are  then  used,  with 
CFL  numbers  shrunk  by  a  factor  (d)-1. 

Systems  of  equations  are  approximated  using  by  now  familiar  field-by-field 
decompositions  ideas.  In  (3.1)  we  replace  the  scalar  constant  a  by  a  constant 
matrix  al,  a  =  {<*..;  where  the  eigenvalues  of  are  non-negative  and  of 

are  non-positive. 

Obviously  a  might  be  taken  to  be  a  sufficiently  large  positive  scalar  multiple  of 
the  identity,  but  this  might  also  lead  to  some  smearing  of  discontinuities  associated 
with  slower  waves.  Other  more  practical  choices  might  involve  freezing  at  some 
constant  state  u,  diagonalizing: 


and  letting 


y~~=T[u) 

ou 


a  =  T{u) 


T-'(u) 


T-l(u) 


where  each  A,  >  |A,(u)|  throughout  the  region.  To  be  safe,  the  margin  of  difference 
between  each  A,  and  the  maximum  value  of  |A,(u)|  has  to  be  sufficiently  large. 
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In  any  case,  the  corresponding  pf(x)  are  obtained  with  the  help  of  the  left 
and  right  eigenvectors  of  §£(«/),  which  we  denote  by  and  i  =  1  We 
interpolate  ■  f±(u)  obtaining  •  pf{x)  exactly  as  in  (3.10)  -  (3.15).  We  then 
define 


.(0 


i-l 


The  fluxes  /*.  are  defined  through  (3.9). 
3 


Generalizations  of  the  type  described  in  (3.20),  (3.21)  using  approximate  Rie- 
mann  solvers  for  j, ,  appropriately  smoothed  at  sonic  points,  may  also  be  ob¬ 
tained. 


Work  is  currently  under  way  with  various  colleagues  applying  these  methods 
to  Euler’s  equations  of  compressible  gas  dynamics  in  multi-space  dimensions. 

IV.  Preliminary  Numerical  Results.  The  numbers  in  this  section  are  often 
written  in  exponential  form,  e.g.  4.2-3  means  4.2  :<  10-3. 


EXAMPLE  1.  The  ENO  schemes  (3.1)  -  (3.15)  in  section  III  combined  with 
the  Runge-Kutta  typeTVD  time  discretizations  (2.19)  -  (2.20)  in  section  II  are  used 
to  solve  the  nonlinear  Burgers  equation  with  periodic  initial  conditions: 


(4-1) 


u(z,0)  =  J  +  5  sin  nx 


-  1  <  x  <  1 


The  exact  solution  is  smooth  up  to  t  =  then  it  develops  a  moving  shock 
which  interacts  with  the  rarefaction  waves.  We  get  the  exact  solution  by  using  a 
Newton  iteration.  For  details,  see  [6]. 


Since  there  is  a  sonic  point,  we  use  the  smooth  LF  (Lax- Friedrichs)  building 
block  in  our  ENO  schemes.  Both  3-3-LF-ENO  (third  order  in  time  and  space  ENO 


schemes  with  Lax-Friedrichs  building  blocks)  and  4-4-LF-ENO  are  used. 


We  use  a  CFL  number  of  0.8  for  3-3-LF-ENO  and  0.6  for  4-4-LF-ENO. 


The  errors  of  the  numerical  solutions  at  t  =  0.3  are  listed  in  Table  1.  Since  the 
exact  solution  is  still  smooth,  we  get  the  full  order  of  accuracy  in  both  Lx  and 
norms. 


At  t  =  the  shock  begins  to  form.  We  use  Ax  =  and  print  out  the  errors 
at  10  points  near  the  shock: 

3- 3-LF-ENO:  -8.7-4,  -2.6-3,  -7.1-3,  -1.3-2,  -1.2-1,  * 

8.2-2,  8.0-3,  1.1-3,  2.6-4,  1.9-4 

4- 4-LF-ENO:  -1.5-4,  -3.2-4,  -1.6-3,  -1.5-2,  -1.2-1,  * 

7.8-2,  7.9-3,  9.2-4,  1.8-4,  9.5-5 


where  the  *  is  the  position  of  the  shock. 


We  see  that  there  is  a  very  good  shock  transition.  (No  oscillations  are  observed). 
Figures  1-6  show  the  shock  transitions. 


In  smooth  regions  the  numerical  solutions  are  very  accurate.  We  compute  the 
Lx  and  L norms  in  the  region  a  distance  of  0.1  from  the  shock  (i.e.  |x—  shock  location|  > 
0.1)  and  list  them  in  Table  2.  From  the  table  we  can  see  that  the  errors  are  of  the 
same  magnitude  as  in  the  smooth  case  when  t  =  0.3. 


At  t  =  1.1,  the  reaction  between  the  shock  and  the  rarefaction  waves  is  over. 
The  solution  becomes  monotone  between  the  shocks.  We  again  print  out  the  errors 
at  10  points  near  the  shock  for  Ax  = 

3-3-LF-ENO:  -1.0-4,  4.6-4,  4.2-4,  -3.3-2.  -8.3-3.  * 
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3.9-2,  1.7-3,  -1.6-4,  -2,5-5,  -7.0-6. 


4-4-LF-ENO:  1.1-5,  6.6-4,  -1.6-3,  -5.7-2,  -1.3-3,  * 
5.7-2,  2.7-3,  -2.4-4,  -5.8-5,  -1.2-6 


Figures  7-12  show  the  shock  transitions. 


The  errors  where  the  solution  is  smooth  are  again  listed  in  Table  2. 


We  can  see  the  excellent  behavior  of  ENO  schemes  in  this  example. 


EXAMPLE  2.  A  two-dimensional  version  of  Example  1 


n +«?).+ of: 

1  u(x,  y,  0)  =  {  +  | 


-  Asinjr(^) 


-  2  <  x,  y<  2 


is  tested  using  the  same  schemes  as  in  Example  1  in  a  dimension  by  dimension 
fashion  (i.e.  T{u)  —  (I  +  A tLx  +  AtLy)(u)  in  (2.1),  together  with  the  R  -  K  time 
discretization:  (2.19)-(2.20)).  The  exact  solution  is  one-dimensional  depending  only 
on  £  =  x  +  y,  however  our  grid  points  are  rectangular  in  (x,  y)  coordinates,  and 
thus  this  example  is  a  truly  2-dimensional  test  problem. 


The  CFL  number  is  always  taken  to  be  half  of  the  one  dimensional  analog,  i.e. 
0.4  for  the  3-3-LF-ENO  and  0.3  for  the  4-4-LF-ENO. 


As  in  Example  1,  we  collect  the  L\  and  Loo  errors  at  t  =  0.3  (smooth  solution) 
in  Table  3  and  the  L i  and  L0 0  errors  in  regions  at  a  distance  of  0.1  from  the  shock 
at  times  t  —  £  and  /  =  1.1  in  Table  4.  We  also  print  out  10  points  near  the  shock 
when  i  =  0,  t  ~  \  and  t  —  1.1  for  Ax  =  Ay  = 
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t  =  J,  x  =  0: 

3-3-LF-ENO:  -9.7-4,  -2.3-3,  -7  6-3,  -4.5-3,  -1.2-1,  * 


8.2-2,  7.9-3,  1.1-3,  2.6-4,  1.9-4 


4-4-LF-ENO:  -1.5-4,  -3.2-4,  -1.6-3,  -1.5-2,  -1.2-1,  * 

7.8- 2,  7.9-3,  9.2-4,  1.8-4,  9.5-5 

t  =  1.1,  x  =  0: 

3- 3-LF-ENO:  -1.0-4,  4.6-4,  4.2-4,  -3.3-2,  -8.3-3,  * 

3.9- 2,  1.7-3,  -1.6-4,  -2.5-5,  -7.0-6 

4- 4-LF-ENO:  1.1-5,  6.6-4,  -1.6-3,  -5.7-2,  -1.3-3,  * 

5.7-2,  2.7-3,  -2.4-4,  -5.8-5,  -1.2-6 

The  shock  transition  graphs  are  very  similar  to  Figures  1-12,  hence  we  omit 
them. 


We  observe  essentially  the  same  results  as  in  the  1-dim  Example  I.  This  indi¬ 
cates  that  our  ENO  schemes  work  well  in  multi-dimensional  problems. 

EXAMPLE  3.  We  use  the  same  schemes  as  in  Example  2  above  to  solve  a 
linear  problem 


(4.3) 


I 


u*  +  ux  +  tiy  =  0,  - 1  <  X,  y  <  1 


u(x<  0) 


1.  H  S 
0,  if  (x,y)£S 


where  S  =  {(x,  y)  :  |x  -  y|  <  |x  +  y|  <  }  is  a  unit  square  centered  at  the 

origin  and  rotated  by  an  angle  of  *  (see  [4]).  We  use  Ax  =  ^  and  run  the  scheme 
up  to  t  =  16  (8  periods  in  time),  in  order  to  study  the  stability  and  the  amount  of 
smearing  of  discontinuities  of  these  methods. 


The  numerical  solutions  at  t  =  2  (after  1  period  in  time)  and  at  t  =  16;  y  =  0 
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and  y  =  —0.4,  are  displayed  in  Figures  13-20. 


Observations :  (1)  the  2  dimensional  schemes  are  stable  under  CFL  numbers 
one  half  of  those  used  for  I  dimension;  (further  experiments  using  1-dimensional 
CFL  numbers  led  to  instability-overflow), 

(2)  The  4-th  order  scheme  resolves  the  discontinuities  better  than  the  3rd  order 
method. 

(3)  Overshoots  and  undershoots,  if  any,  are  negligible. 

EXAMPLE  4.  We  did  not  prove  for  these  ENO  methods  that  limit  solutions 
satisfy  the  entropy  condition.  However,  numerical  experiments  in  [6],  including 
some  tests  using  nonconvex  fluxes,  indicated  the  convergence  of  ENO  schemes  to 
the  correct  entropy  solution.  We  test  our  schemes  3-3-LF-ENO  and  4-4-LF-ENO 
for  Riemann  problems  for  two  such  fluxes.  One  is 

((4.4))  /<«)  =  i(«3  -  l)(«l  -  4) 

4 

with  «£,  =  2,  Ur  =  —2  (the  exact  solution  is  a  shock  followed  by  a  rarefaction 
wave  followed  by  another  shock)  and  with  «£,  =  —3,  ur  =  3  (a  stationary  shock  at 
x  =  0).  See  [6)  for  details.  Our  schemes  converge  to  the  correct  solutions  in  both 
cases  with  good  resolution.  The  results  are  displayed  in  Figures  21-32. 


Another  nonconvex  flux  we  test  is  the  well  known  Buckley-Leverett  example: 


/(«)  = 


4u3  +  ( 1  -  u)3 


with  initial  data  u  =  1  in  [  —  £,0],  and  u  =  0  elsewhere. 


The  exact  solution  is  a  shock-rarefaction-contact  discontinuity  mixture.  Our 
schemes  resolve  the  correct  solution  well.  However,  the  3-3-LF-ENO  (using  the  Lax- 
Friedrichs  building  block)  smears  more  than  the  3-3-EO-ENO  (using  the  Engquist- 
Osher  building  block)  (Figures  33-38).  Although  in  this  example  /'(«)  >  0  so 


there  was  no  need  to  smooth  the  flux,  the  observed  improvement  using  the  upwind 
building  block  indicate  that  sometimes  it  is  worthwhile  spending  the  effort  to  smooth 
the  EO  or  other  upwind  flux  rather  than  to  use  the  simple  LF  building  block. 
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Table  3.  (Example  2). 

t  =  0.3;  E  :  type  of  error;  r  :  numerical  order 
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Figure  36:  3-3-EO-ENO,  Ax 
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